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IMPLEMENTATION  OF  THE 

LAMONT-DOHERTY  GEOLOGICAL  OBSERVATORY  NORMAL  MODE/FAST  FIELD  MODEL 
ON  THE  NUSC  VAX  780/11  COMPUTER 


1.  INTRODUCTION 


A  common  boundary  value  problem  in  underwater  acoustics  is  to  determine 
the  sound  pressure  field  caused  by  a  harmonic  point  source  located  within  a 
plane  multilayered  medium.  The  classic  solution  to  this  problem,  when  cast  in 
a  cylindrical  coordinate  system,  takes  the  form  of  a  Hankel  transform.1  This 
integral  solution  can  be  manipulated  to  produce  several  different 
representations  for  the  field2  from  which  numerical  evaluations  can  be 
obtained.3  A  recent  survey4’5  lists  existing  computer  models  based  on  these 
representations  and  aimed  at  obtaining  predictions  of  sound  propagation  for 
use  in  sonar  applications. 

Most  of  the  models  in  the  above  survey  deal  with  the  propagation  of 
pressure  waves  within  layered  liquid  media.  An  exception  is  found  in  the 
models  developed  by  Kutschale,  which  account  for  both  pressure-wave  and 
shear-wave  propagation  within  layered  liquid/solid  media.6’7  Kutschale's 
models  were  originally  developed  for  application  in  the  Arctic  Ocean,  where 
the  combination  of  a  solid  ice  canopy  and  an  upward  refracting  ocean  sound 
speed  profile  requires  the  capability  for  modeling  the  interaction  between  the 
two  types  of  waves.  More  recently,  for  shallow  water  environments,  there  has 
been  interest  in  modeling  the  loss  of  acoustic  energy  into  the  ocean  bottom  as 
a  result  of  pressure-wave  to  shear-wave  conversion.  Of  course,  geophysicists 
have  always  been  concerned  with  the  propagation  of  seismic  waves  in  layered 
sol  id  media.9 ’10 

This  report  describes  the  implementation  of  Kutschale's  LDGO  normal 
mode/fast  field  model  on  the  NUSC  VAX  780/11  computer.  The  version  obtained 
by  NUSC,  designated  17HH,  permits  two  computational  methods  for  evaluating  the 
field  due  to  a  point  source  located  within  a  layered  fluid/solid  medium:  (1) 
a  normal  mode  method,  in  which  the  field  is  represented  as  a  finite  sum  of 
propagating  modes  (discrete  spectrum)  plus  a  branch  line  integral  (continuous 
spectrum),  and  (2)  a  fast  field  method,  in  which  the  field  is  represented  in 
the  form  of  a  finite  Fourier  integral  whose  evaluation  can  proceed  directly 
via  a  fast  Fourier  transform  (FFT).  In  addition  to  the  capability  for 
treating  propagation  of  both  compressional  and  shear  waves,  the  17HH  model  can 
accommodate  the  effects  of  rough  boundaries  via  new  boundary  conditions 
derived  from  perturbation  analysis.11’12  We  designate  the  NUSC  version  of  the 
computer  code  as  NUSC17HH. 

Section  2  summarizes  the  basic  theory  leading  to  the  two  representations 
of  the  field  used  by  NUSC17HH  and  presents  some  of  the  computational 
limitations  of  the  normal  mode  option  that  must  be  considered  when  running  the 
model.  Section  3  presents  a  detailed  description  of  the  data  input  sequence 
required  to  obtain  NUSC17HH  predictions  of  sound  propagation  and  the  data 
input  sequence  that  must  be  supplied  to  a  locally  developed  plotting  program 
(PL0T17HH).  Section  4  describes  a  sample  command  file  structure  needed  to  run 
the  model  as  a  batch  job  submission.  Finally,  in  Section  5  we  present  several 
numerical  examples  that  illustrate  various  features  and  capabilities  of  the 
model . 


1 


TR  7221 


2.  BASIC  THEORY 


2.1  MATHEMATICAL  MODEL 

Let  the  half-space  z  >  0  (z  positive  down)  of  a  cylindrical  coordinate 
system  (r,  9,  z)  be  occupied  by  a  layered  elastic  medium  as  shown  in  figure  1. 
Each  layer  is  characterized  by  five  material  properties,  i.e.,  compressional 
wave  sound  speed  (c),  shear  wave  sound  speed  (v),  compressional  wave 
attenuation  (a),  shear  wave  attenuation  (3),  and  density  (p).  The  parameters 

v.,  a.,  3  ,  and  p.  are  presumed  to  be  constant  within  the  j*"  layer.  The 
J  J  J  J 

compressional  sound  speed,  c.,  is  presumed  to  vary  with  depth  according  to 

J 

l/c?(z)=a .+b.z.  The  layer  constants  a.  and  b.  are  determined  by  the  values  of 

J  J  J  J  J 

Cj  at  Zj+0  and  Zj+^-0.  The  region  z  >  z^  is  taken  to  be  a  homogeneous  solid 

(or  liquid)  whereas  the  region  z  <  0  is  taken  to  be  a  homogeneous  fluid  (or 
vacuum).  In  this  model,  layer  interfaces  are  presumed  to  occur  at  both  the 
source  and  receiver  depths.  Any  discontinuities  in  material  properties  of  the 
medium  can  be  accommodated  at  a  layer  interface. 

2.2  BASIC  EQUATIONS 

With  the  passage  of  a  harmonic  (e“1a)t,  uj=2tt f )  elastic  wave  through  the 

JL,  U 

jTn  layer,  the  resulting  particle  displacement  within  z.  <  z  <  z.  ,  can  be 

J  J  ^ 

found  by  solving  the  reduced  wave  equations  satisfied  by  the  two  displacement 
potentials,  <J>.  and  4'i,9  i.e., 

J  J 


(1) 


and 


(2) 


where 


V2  =  Ii_(r  i_)  +li 


r  3r  3r  3z2 


(3) 


The  vertical  (q.)  and  horizontal  (w-)  components  of  displacement  are  recovered 

J  vJ 


from  these  potentials  by  differentiation,  i.e., 


(4) 


2 


0 


layer  0 

layer  1 

• 

• 

• 

layer  j 

• 

© 

• 

layer  N 

1 

1  z  layer  N+l 

Z1 = 
z2 


z 


j 


z 


j+1 


ZN 

ZN+1 


v(z)  =  v(Zj+0) 
p(z)  =  P (Zj+0) 
a(z)  =  a(zj+0) 
3(z)  =  e(zj+0) 


/ 

/ 

/  c“2(z)  =  a.  +  b  z 

J  J  J 

/ 

/ 

/ 

Zj+1 


Figure  1.  The  Layered  Model  Showing  (a)  Layer  Index  Notation, 
and  (b)  Variation  of  Acoustic  Properties  within  Zj  <  z  <  Zj+^ 
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and 


3<j>  .  32\|j  .  go2 

W.  =  — ]  + - J  +  —♦_••  (5) 

3  3  z  3  z2  v2.  3 

j 

Similarly,  the  normal  (p  )  and  tangential  (p  )  components  of  stress  are 
obtained  via 


and 


zz 


u/  3  w . 

X  .  —  <|> .  +  2p  .  — -J  , 


3  c2  J 


3  Z 


rz 


3q.  3  w • 

-^(-J  +-J)  . 


3  z  3  r 


(6) 


(7) 


Here  X.  and  y  .  are  Lame's  elastic  parameters  which  are  related  to  the  sound 

J  J 

speeds  and  density  by 


and 


=  xj 


+  2y 


(8) 

(9) 


Equations  (1)  and  (2)  are  conveniently  solved  using  Hankel  transforms. 
The  displacement  potentials  <t>  i  and  'I',  can  be  represented  in  the  form 

J  J 


<t>  j (r,z)  =  /  fj(k,z)J0  (kr)kdk  ,  (10) 

o 

and 

oo 

tj(r,z)  =/  gj(k,z)J0  (kr)kdk  ,  (11) 

o 

where  Jq  is  the  Bessel  function  of  the  first  kind,  k  is  the  horizontal 
wavenumber,  and  f  and  g.  are  solutions  to  the  depth-dependent  equations 

J  J 


4 
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d2  f .  to2 

— J  +  (—  -  k2)fi  =  0  ,  (12) 

dz^  <£  J 

sJ 

and 

d2  g .  u)2 

— J  +  (—  -  k2)9i  =  0  .  (13) 

dz2  v?  J 


If  c.  and  v,  are  constant  within  a  layer,  then  both  f  and  g  are  linear 
J  J  J  J 

combinations  of  complex  exponential  functions  representing  upgoing  and 
downgoing  waves.  When  1/c^  is  a  linear  function  of  z,  then  the  appropriate 

solution  to  equation  (12)  comprises  a  linear  combination  of  Airy 
functions.1  ,3 ,15  Each  f,  and  g  involves  two  arbitrary  constants  that  must  be 

J  J 

determined  for  each  layer.  For  j  =  0  and  j  =  N+l,  the  potentials  <j>  .  and  <|> . 

J  <1 

must  tend  to  zero  as  |z|  -*■  °°.  The  remaining  conditions  for  determining  the 
constants  are  provided  by  requiring  the  solutions  to  equations  (12)  and  (13) 
in  each  layer  to  satisfy  appropriate  boundary  conditions  at  layer  interfaces. 
For  elastic  media  in  "welded  contact,"  the  components  of  displacement,  q.  and 

si 

w.,  and  stress,  (p  ).  and  (p  ).,  remain  continuous  across  the  planes  z  =  z. 

3  zzjrzj  j 

and  z  =  zj+i*  These  continuity  requirements  are  sufficient  to  determine  the 

unique  solutions  to  the  depth-dependent  equations  for  layered  media.  Although 
equations  (1),  (2),  (12),  and  (13)  hold  only  in  source-free  regions,  the 
inclusion  of  a  point  acoustic  source  can  be  accommodated  by  requiring  that  the 
displacement  potential  s  <t>  •  andti  satisfy  additional  boundary  conditions  on 

the  plane  z  =  zs.9 

For  z  >  zs,  each  displacement  potential  for  the  layered  medium  can  be 
cast  in  the  canonical  form16  ,17 


<}>(r,z)  =  /  rfir}  U(zs,k)V(z,k)J0  (kr)kdk  ,  (14) 

o 

and 


*(r,z)  =  /  U(zs,k)V(z,k)J0(kr)kdk  ,  (15) 

o 


where  the  multipliers  of  J0(kr)kdk  represent  the  Green's  function  solutions  to 
the  systems  of  depth-dependent  equations  and  can  be  determined  from  the  layer 
solutions  f.  and  g.,  j  =  0,  1,  ...,  N+l  using  the  classic  Thomson-Haskel 1 

kJ  J 


5 
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matrix  approach.5’6’15  Here  U(z,k)  represents  the  solution  to  equation  (12) 
that  satisfies  all  boundary  conditions  for  z  <  z$,  and  V(z,k)  represents  the 

solution  that  satisfies  all  boundary  conditions  for  z  >  z.^.  When  z  <  z$,  the 

appropriate  form  of  each  potential  is  obtained  by  interchanging  z  and  z$  in 

equations  (14)  and  (15).  The  denominator  A (k)  =  p(z$)W(k)  where  W(k)  is  the 

Wronskian  of  U  and  V,  i.e.,  W(k)  =  V(zg  >k)Uz(zs  ,k)  -  Vz (z$ ,k)U(zs  ,k) .  These 

Green's  functions  account  for  all  the  environmental  properties  of  the  layered 
medium,  as  well  as  the  locations  of  the  source  and  receiver.  Their  evaluation 
for  any  problem  constitutes  the  major  computational  effort  required  to  obtain 
predictions  of  sound  propagation. 

2.3  NORMAL  MODE  REPRESENTATION 

Integrals  of  the  form  in  equations  (14)  and  (15)  are  improper.  For 
lossless  media,  the  Green's  functions  appearing  in  the  integrands  are  singular 
at  (real)  roots  of  A(k).  Using  the  contour  integral  techniques  presented  by 
Ewing,  Jardetzky  and  Press,9  the  modal  representation  for  the  displacement 
potential  <j>  is  given  by 


M  ,  v 

<t>(r,z)  =  l  C0(km)U(zs,km)U(z,km)H011  J(kmr) 
m=l 

+  /  Ci  (k)U(z  ,k)U(z,k)H0  ^  )(kr)kdk  ,  (16) 

BL 

with  a  similar  representation  for  4).  The  discrete  wavenumbers  km,  m  =  1,  ..., 
M  are  roots  of  the  period  equation,  A(k)=0.  For  lossless  liquid  media,  the 


Im(k) 


continuous 
spectrum  «■ 
region 


xxxxxxxxxxxxxxxx 


discrete 
«•  spectrum 
region 


exponential ly 
decaying 
region 


+ 

k=o>/c 


poles 


N+l 


k=o)/min(c, ,  j  =  l,. . .  ,N) 

J 


Re(k) 


Figure  2.  Sketch  of  the  Complex  Wavenumber  Plane  Showing  the  Location  of  the 
Discrete  Mode  Eigenvalues  (poles),  the  Branch  Cut  (xxxxx)  and  the 
Contour  ( .  )  for  the  Branch  Line  Integral 
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location  of  these  poles  and  the  branch  line  contour  (BL)  is  depicted  in  figure 
2.  For  elastic  media,  two  branch  line  contours  are  obtained,  associated  with 
the  branch  points  at  k=^/c^+j  and  k^/v^.  Each  term  appearing  in  the  sum  is 

termed  a  normal  mode.  The  factor  U(zs>km)  determines  the  influence  of  the 

source  on  each  mode  amplitude.  The  vertical  distribution  of  each  mode  is 
described  by  the  factor  U(z,km).  The  coefficient  Co  gives  the  relative 

excitation  of  each  mode.  In  the  case  of  liquid  layers  for  which  l/c£  varies 

J 

linearly  with  depth  z,  explicit  expressions  for  the  coefficients  Co  and  Ci  may 
be  found  in  Stickler's  paper.15 

2.4  FAST  FIELD  REPRESENTATION 

Although  equations  (14)  and  (15)  are  improper,  the  significant 
contribution  to  these  integrals  occurs  within  a  finite  range  of  the  wavenumber 
axis,  ko  <  k  <  kmax.  This  suggests  that  these  integrals  may  be  evaluated  by 

numerical  quadrature.  An  efficient  method  for  carrying  out  the  numerical 
integrations  is  possible  if  the  integrals  can  be  manipulated  into  the  form  of 
finite  Fourier  integrals  so  that  the  FFT  algorithm  can  be  used.  In  underwater 
acoustics  applications,  this  approach  is  known  as  the  fast  field  method.1’3’7 

To  this  end,  we  use  the  identity  20o  (kr)  =  Ho^^(kr)  -  Hg^^(-kr)  to  cast 
equation  (14)  in  the  equivalent  form 


<t>  (r ,z)  =  J  D(z,zs;k)Ho  h  )(kr)kdk. 


(17) 


where  we  have  set  D(z,zs;k)  =  A(k)[2A (k)]"1  U(z  ,k)U(z,k).  Next,  for  kr  »  1, 

we  replace  the  Hankel  function  Ho^^(kr)  by  the  first  term  in  its  asymptotic 
expansion,  i.e.. 


(18) 


Under  these  conditions,  the  integral  appearing  in  equation  (17)  may  be 
approximated  by 


k 

max 

4>(r,z)  =  (2/irir)1/2  /  D(z,z  ;k)  exp(ikr)  k1/2  dk  , 

ko 


N-l  km+l 

=  (2/Trir)172  l  j  D(z,z  ;k)  exp(ikr)  k1/2  dk,  (19) 

m=0  k 

m 
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where  the  discrete  wavenumber  variable  k  =  kn  +  mAk,  m=0,  1,  ....  N-l  has 

m  u 

been  introduced  spanning  the  range  kmax  -  k0  =  (N-l)Ak.  Provided  the  sampling 

interval  Ak  is  sufficiently  small,  the  integrals  appearing  in  equation  (19) 
may  be  evaluated  approximately  via  the  Trapezoid  Rule  leading  to  the  result 


N-l 

4>(r,z)  *  Ak(2/iTir ) 1/2  l  D(z,z  ;k  )  exp(ikmr)  ki/2  .  (20) 

m=0 


If  the  range  interval  rmax  -  Tq  =  (N-l)Ar  is  similarly  sampled  according  to  r 

=  ro  +  nAr,  n=0,  1,  ...,  N-l  subject  to  the  FFT  constraint  ArAk  =  2tt/N,  the 
displacement  potential  of  equation  (20)  has  the  following  finite  complex 
sequence  approximation: 


N-l 

<t>(rn,z)  =  Ak ( 2/TTi  rn ) 1 72  exp(ik0rn)  l  Em  exp(2itimn/N) .  (21) 


The  above  series  has  the  desired  form  of  a  discrete  Fourier  series  that  can  be 
evaluated  rapidly  using  the  FFT  algorithm  provided  that  N  is  a  power  of  2. 

The  input  sequence  Em  is  determined  to  be 

Em  =  km/2  D^z>zs;km)  exP(imAkro).  (22) 

The  result  of  letting  the  complex  sequence  Em,  m=0,  1,  ...,  N-l  be  input 

to  a  complex  FFT  routine  is  the  complex  sequence  <J>(r  ,z)  at  N  distinct  range 

points  rn,  n=0,  1,  ...,  N-l.  Only  the  first  N/2  values  are  considered  valid 

due  to  aliasing.  It  follows  that  N  estimates  of  Em  and  consequently 

D ( z , zs ; km )  as  a  function  of  the  sampled  wavenumber  km  are  required.  The  time 

consuming  portion  of  the  fast  field  method  lies  in  the  generation  of  this 
sampled  depth-dependent  Green's  function. 

2.5  DISCUSSION 

The  normal  mode  and  fast  field  computational  options  for  NUSC17HH  are 
summarized  in  figure  3.  The  program  flow  for  these  choices  is  directed  by  the 
two  flags  IM0DE  and  IC0NT.  The  value  IC0NT  =  0  limits  the  normal  mode 
representation  to  the  discrete  mode  sum  only.  Although  the  diagram  suggests 
complete  symmetry  for  the  two  methods  of  representing  the  acoustic  field,  the 
normal  mode  implementation  is  limited  by  the  following  constraints: 

1.  Only  flat  boundaries  can  be  treated.  The  rough  boundary  conditions 
developed  by  Kuperman  and  Ingenito11 512  have  not  been  incorporated  into  the 
modal  representation. 
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2.  Only  lossless  layers  can  be  treated.  This  implies  that  and  3^ 
must  be  set  to  zero  for  all  j,  j=0,  . ..,  N+l  in  the  normal  mode  method. 

3.  Only  the  branch  line  integral  associated  with  the  compressional  sound 

speed  cN+1  is  included  if  the  flag  IC0NT  *  0  is  selected.  Moreover,  only  that 

portion  of  the  branch  line  integral  along  the  real  axis  is  computed.  The 
contribution  of  the  branch  line  integral  associated  with  the  shear-wave  sound 
speed  of  the  half-space  z  >  is  not  evaluated  in  the  17HH  version  of  the 

model.  In  addition,  in  the  evaluation  of  the  branch  line  integral,  the 
variation  of  sound  speed  within  each  layer  is  taken  to  be  constant. 

4.  If  IC0NT  <  0,  the  sampling  interval  Ar  and  maximum  range  NAr/2  which 
is  automatically  determined  for  the  branch  line  contribution  is  not  usually 
commensurate  with  the  values  of  DELR  and  RMAX  selected  by  the  user  (see 
NAMELIST  &CARD4  of  the  next  Section).  In  this  case,  the  sampling  interval  is 
determined  by  the  distribution  of  N  equispaced  values  of  k  around  the  branch 
cut  on  the  real  axis. 

The  fast  field  method  of  NUSC17HH  offers  the  prospective  user  the  most 
general  option  for  obtaimnq  numerical  predictions  of  sound  propagation  within 
a  layered  liquid/solid  medium.  With  this  choice  of  computational  method,  both 
rough  boundaries  and  absorbing  media  can  be  accommodated.  In  addition,  the 
effects  of  the  continuous  spectra  associated  with  both  the  compressional  waves 
and  shear  waves  can  be  evaluated. 
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Wave  Equation 
+ 

Boundary  Conditions 


Integral  Solution  for  a 
Harmonic  Point  Source  in  a 
Layered  Fluid/Solid  Media 
(Hankel  Transform) 


Normal  Mode  Method 


Fast  Field  Method 


( IM0DE  =  1) 


(IM0DE  =  0) 


Contour  Integral 
Representation 


Fourier  Integral 
Representation 


Proper  Mode  Sum  + 
(Discrete  Spectrum) 


Branch  Line  Integral 
(Continuous  Spectrum) 


(IC0NT  =  0) 


( IC0NT  *  0) 


Direct  Numerical 
Evaluation  via  FFT 


Fast  Field  Method 
Eval uation 


Gaussian  Quadrature 
Eval uation 


( IC0NT  <  0) 


( IC0NT  >  0) 


Figure  3.  Basic  Computational  Options  for  NUSC17HH 
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3.  DATA  INPUT  SEQUENCES 


3.1  DATA  SEQUENCE  FOR  NUSC17HH 

The  data  input  sequence  for  NUSC17HH  has  been  restructured  from  that  of 
the  LDGO  version  to  make  use  of  the  FORTRAN  NAMELIST  feature.  In  addition, 
the  input  variables  have  been  renamed  to  make  their  function  more  apparent  to 
the  casual  user.  An  equivalence  list,  which  relates  the  NUSC  input  sequence 
names  to  the  LDGO  sequence  names,  is  provided  in  the  appendix. 


CARD  1  —  NAMELIST  &CARD1  inputs  N0SSP,  IM0DE,  ITEST,  IPL0T 


Variable  name  Description 


N0SSP  No.  of  points  In  the  sound  speed  profile. 

Each  point  corresponds  to  a  depth  z  at  which 
the  val ues  of  c,  v,  a,  B,  p  and  o  are 
supplied.  Values  at  both  the  source  and  the 
receiver  depths  must  be  included  in  the 
profile.  See  CARD  2. 

IM0DE  Flag  to  select  the  method  of  computation. 

IM0DE  =  0  gives  fast  field  method 
=  1  gives  normal  mode  method 

ITEST  Flag  to  select  the  internal  test  case 

profile.  All  default  parameters  for  the 
test  case  correspond  to  those  used  in  the 
AESD  Workshop  Case  IB.16 

ITEST  =  0  gives  user-supplied  profile 
=  1  gives  internal  profile 

IPL0T  Flag  to  save  the  kernel  and  propagation  loss 

values  (fast  field  output)  or  excitation 
function  and  propagation  loss  values  (normal 
mode  output)  in  the  file  17HH.PLT  for 
subsequent  access  by  the  plotting  program 
PL0T17HH. 

IPL0T  =  0  gives  no  plot  file 
=  1  gives  plot  file 
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CARD  2  —  inputs  ((V(I,J),  J=1 ,7 ) , 1=1 ,N0SSP)  -  FORMAT  (7D10.0) 


Variable  Name  Description 


V(I,J)  Array  of  physical  properties  of  the 

environmental  profile.  The  values  at  the 
source  and  receiver  depths  must  be  included 
in  the  profile.  Each  point  in  the  profile 
separates  layers  having  different  acoustic 
properties.  Discontinuities  in  acoustic 
properties  can  be  accommodated  by  assigning 
two  sets  of  values  at  the  same  depth.  For 
each  I-index  profile  point,  the  7  physical 
values  needed  are 


V(I.l) 

z 

(■). 

Depth 

V( I ,2) 

c 

(m/s) , 

Pressure-wave  sound  speed 

V( I ,3) 

V 

(m/s) , 

Shear-wave  sound  speed 

V(I»4) 

a 

(dB/km) , 

Pressure-wave  attenuation 

V  ( 1 , 5 ) 

3 

(dB/km) , 

Shear-wave  attenuation 

V(I,6) 

P 

(g/cm3) , 

Density 

V  ( I » 7 ) 

a 

(■). 

Rms  roughness 

Note:  In  the  LDGO  version  of  the  code,  the 
normal  mode  computational  method  was 
selected  if  V( 1 ,4)  =  V(  1 ,7 )  =  0  for  1  <  I  < 
N0SSP,  and  if  REFLTO  =  0  (see  CARD  5).  In 
the  NUSC  version,  if  IM0DE  =  0,  and  the 
profile  attenuations  and  surface  roughness 
are  everywhere  zero,  then  V( 1 ,4 )  is  set 
equal  to  Ak  at  all  depths.  (Ak  is  defined  in 
RMIN  in  CARD  4.)  This  attenuation  is 
subsequently  removed  from  the  solution. 
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Variable 

ISRCE 

IRCVR 

IPRUN 


N0FFT 

INEXT 


ITAPR 

IDISP 


IC0NT 


Name  Description 


Profile  I-index  of  the  source. 

Profile  I-index  of  the  receiver. 

Pruning  factor  to  limit  display  of  aliased 
propagation  loss  values.  Only  the  first 
1/ IPRUN  of  the  propagation  loss  versus  range 
curve  is  displayed.  This  parameter  applies 
to  the  fast  field  method  only.  IPRUN  >  2. 

No.  of  points  of  the  sampled  integrand  to  be 
used  as  input  to  the  inverse  FFT.  This 
parameter  applies  to  the  fast  field  method 
only.  N0FFT  <  8192. 

Flag  to  select  the  entry  point  for  the  next 
data  input  sequence. 

INEXT  =  0  gives  an  exit  after  this  case 

<  0  gives  an  entry  point  at  CARD2 

>  0  gives  an  entry  point  at  CARD3 

Flag  to  select  the  option  of  tapering  the 
integrand  before  taking  the  inverse  FFT. 

ITAPR  =  0  gives  no  tapering  (square  window) 
=  1  gives  tapering  at  both  ends 
Flag  to  select  the  field  variable  of 
interest  for  this  data  set. 

IDISP  =  0  gives  pressure  as  output 

<  0  gives  horizontal  displacement 

>  0  gives  vertical  displacement 

Flag  to  select  the  computational  method  for 
evaluating  the  branch  line  integral.  This 
applies  to  the  normal  mode  method  only. 

IC0NT  =  0  gives  discrete  modes  only 

<  0  gives  fast  field  evaluation 

>  0  gives  Gaussian  quadrature 

eval uation 
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CARD  4  —  NAMELIST  &CARD4  inputs 

FREQ,  SLVL,  CMIN,  CMAX, 

RMIN,  RMAX,  DELR 

Variable  Name 

Description 

FREQ 

Source  frequency  (Hz) 

SLVL 

Source  level  (dB  re  1  uPa  @  1  m) 

CMIN 

Sound  speed  minimum  (m/s).  id/CMIN  gives  the 
highest  wavenumber  in  the  integration  for 
the  fast  field  method.  For  the  normal  mode 
method,  it  supplies  the  upper  wavenumber 
limit  for  the  root  finding  routine  for  the 
mode  eigenvalues.  If  CMIN  =  0,  a  default 
value  will  be  determined  from  the  sound 
speed  profile  for  the  normal  mode  method. 

CMAX 

Sound  speed  maximum  (m/s).  w/CMAX  gives  the 
lowest  wavenumber  in  the  integration  for  the 
fast  field  method.  For  the  normal  mode 
method,  it  supplies  the  lower  wavenumber 
limit  for  the  root  finding  routine  for  the 
mode  eigenvalues.  If  CMAX  =  0,  a  default 
value  will  be  determined  from  the  sound 
speed  profile  for  the  normal  mode  method. 

In  addition,  the  integration  will  start  from 
zero  for  the  fast  field  method. 

RMIN 

The  starting  range  (in  km)  of  the 
propagation  loss  curve  for  the  normal  mode 
method.  For  the  fast  field  method,  this 
range  (in  m)  is  given  by 

Ar  =  2ir/(N0FFT*Ak) 
where 

(NQFFT-l)Ak  =  ( u/CMIN  -  o>/CMAX). 

RMAX 

The  maximum  range  (in  km)  of  the  propagation 
loss  curve  for  the  normal  mode  method.  For 
the  fast  field  method,  this  range  (in  m)  is 
given  by  N0FFT*Ar. 

DELR 

The  range  step  (in  km)  of  the  propagation 
loss  curve  for  the  normal  mode  method.  For 
the  fast  field  method,  this  step  (in  m)  is 
given  by  Ar  =  RMIN. 
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CARD  5  --  NAMELIST  &CARD5  inputs 

REFLTO,  ALPHAO,  ATTENO,  RH0O,  H 

Variable  Name 

Description 

REFLTO 

Rms  surface  roughness  (m)  at  z  =  0  . 

If  REFLTO  *  0,  then  the  next  three 
parameters  are  internally  set  to  zero. 

ALPHAO 

Pressure-wave  sound  speed  (m/s)  of  the  fluid 
in  the  upper  half-space  z  <  0.  This 
parameter  determines  whether  this  region  is 
a  vacuum  or  a  fluid. 

ATTENO 

Pressure-wave  attenuation  (dB/km)  of  the 
fluid  in  the  upper  half-space  z  <  0. 

RH0O 

Density  (g/cm3)  of  the  fluid  in  the  upper 
half-space  z  <  0. 

H 

Separation  distance  (m)  between  the  source 
and  receiver  when  located  within  the  fluid 
half-space  z  <  0. 

To  use  this  option,  the  following 
requirements  must  be  met: 

1.  V ( 1 , 2 )  =  V(2,2)  =  0.  See  CARD2. 

2.  V( 1+1,2)  -  V( 1 ,2)  *  0  for  at  least  one 
value  in  the  range  2  <  I  <  N0SSP  -  1. 

3.  V( I ,3)  =  0  for  1  <  I  <  N0SSP. 

4.  Both  ALPHAO  and  H  must  be  specified. 
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3.2  DATA  SEQUENCE  FOR  PL0T17HH 

A  plotting  program  (PL0T17HH)  was  developed  to  provide  quality  plots  of 
the  curves  generated  by  NUSC17HH  and  stored  in  the  output  file  named  17HH.PLT. 
Plotting  is  effected  online  so  that  prompts  for  data  input  are  issued  to  a 
user's  terminal.  Plots  can  be  directed  to  several  display  devices,  including 
the  Tektronix  4014,  Cal  comp  1039,  and  FR80.  An  explanation  of  the  input 
prompts,  which  allow  some  flexibility  in  plotting  the  data,  are  given  below. 


Variable  Name  Description 


IBFIL 

Flag  to  select  plot  device. 

IBFIL  =  2  to  select  Tektronix  4014 
=  4  to  select  FR80 
=  8  to  select  Cal  comp  1039 

ISP0T 

Flag  to  select  line  thickness.  This  flag  is 
used  only  for  plots  directed  to  the  FR80. 

ISP0T  =  0-7  where  the  line  thickness 
increases  with  ISP0T 

IM0DE 

Flag  to  identify  the  method  of  computation 
used  in  NUSC17HH. 

IM0DE  =  0  signifies  the  fast  field  method 
=  1  signifies  the  normal  mode  method 

IUNIT 

Flag  to  select  the  units  of  range  and  dB. 

IUNIT  =  0  gives  km  and  dB  re  1  m 
=  1  gives  kyd  and  dB  re  1  yd 

RMIN 

Minimum  range  displayed  on  the  range  axis. 

RMAX 

Maximum  range  displayed  on  the  range  axis. 

RPIN 

No.  of  range  units  per  division  on  the  range 
axi  s . 

TMIN 

Minimum  dB  level  on  the  propagation  loss 
axis. 

TMAX 

Maximum  dB  level  on  the  propagation  loss 
axis. 

TPIN 

No.  of  dB  per  division  on  the  propagation 
loss  axis. 
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4.  COMMAND  FILE  SEQUENCE 


The  sample  command  file,  designated  NUSC17HH.C0M,  can  be  used  to  place 
the  program  NUSC17HH  into  the  batch  queue.  The  command  file  comprises  the 
following  lines: 


SASSIGN  NL  SYS$PR I  NT 
$ ASS  I GN  17HH.DAT  F0ROO5 
SASSIGN  17HH.PLT  F0ROO7 
$RUN  DBAO : [TH0MS0N]NUSC17HH 
$C0PY  NUSC17HH.L0G  17HH.0UT 
SEXIT 


line  1 
line  2 
1  ine  3 
line  4 
1  ine  5 
line  6 


Line  1  -  Assigns  the  log  file  to  the  user's  disc  area. 

Line  2  -  Assigns  the  data  file  17HH.DAT  to  unit  5. 

Line  3  -  Assigns  the  plot  file  17HH.PLT  to  unit  7 

for  subsequent  access  by  the  program  PL0T17HH. 

Line  4  -  Runs  the  executable  version  of  the  program, 
i.e.,  NUSC17HH.EXE. 

Line  5  -  Copies  the  log  file  which  contains  all  printed 
output  from  NUSC17HH  (written  to  unit  6) 
to  the  output  file  17HH.0UT 
Line  6  -  Terminates  the  command  file. 

The  command  file  is  invoked  by  keying  in  "SUBMIT  NUSC17HH".  The  plotting 

program  PL0T17HH  is  invoked  by  keying  in  "RUN  DBAO:[TH0MS0N]PL0T17HH"  and 

responding  to  the  prompts  described  in  section  3.2. 
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5.  NUMERICAL  EXAMPLES 


In  this  section,  we  present  several  numerical  examples  to  illustrate  the 
computational  capabilities  of  NUSC17HH  and  the  graphical  capabilities  of 
PL0T17HH.  The  printed  output  of  the  model,  which  is  stored  in  the  file 
17HH.0UT,  will  not  be  described.  However,  since  some  of  the  printed  output 
relates  to  the  variable  names  used  by  the  LOGO  version  of  the  model,  the  user 
may  find  it  convenient  to  refer  to  the  equivalence  list  given  in  the  appendix. 

5.1  TEST  CASE  1. 

The  first  example  (test  case  1)  is  taken  from  an  AESD  Workshop  on 
acoustic  modeling.16  The  data  input  sequence  for  this  example  is  shown  in 
figure  4.  By  setting  the  flag  ITEST  =  1,  all  the  parameters  listed  in  figure 
4  are  set  automatically  via  DATA  statements,  so  that  an  equivalent  data  input 
sequence  is  given  by  (note  that  each  appears  in  column  2) 

1  2  3  4  5  6  7 

123456789012345678901234567890123456789012345678901234567890123456789012345678 

&CARD1 

ITEST=1 , 

&END 

&CARD3 

SEND 

&CARD4 

&END 

&CARD5 

&END 

The  canonical  North  Pacific  Ocean  sound  speed  profile  for  this  case  is 
depicted  in  figure  5.  The  ocean  bottom  comprises  a  uniform  half-space  of 
density  1.9176  g/cm3  and  a  sound  speed  of  1555.52  m/s.  Calculations  of  the 
acoustic  field  were  carried  out  at  25  Hz  for  a  source  at  a  depth  of  253.9  m 
and  a  receiver  depth  of  863.5  m. 

Figure  6  presents  the  results  for  test  case  1  using  the  normal  mode 
computational  option.  The  discrete  spectrum  comprised  44  propagating  modes. 
The  contribution  of  the  continuous  spectrum  was  not  included  by  setting  IC0NT 
=  0  and  limiting  the  search  area  of  the  root  finding  routine  to  the  range  of 
wavenumbers  within  w/CMIN  =  oj/  147 1 . 88  <  k  <  cu/ 1555 . 52  =  u/CMAX  determined  by 
the  maximum  and  minimum  values  of  the  sound  speed  in  the  profile.  Figure  6(a) 
shows  the  normalized  magnitude  of  the  excitation  function  of  each  mode  in  the 
sum.  Figure  6(b)  shows  the  propagation  loss  versus  range  curve,  a  coherent 
sum  of  44  modes,  which  exhibits  a  well  defined  convergence  zone  behaviour  for 
this  example. 

Figure  7  shows  the  results  obtained  for  test  case  1  using  the  fast  field 
computational  option.  For  these  results,  512  points  were  used  to  sample  the 
discrete  mode  region  of  the  wavenumber  axis.  Figure  7(a)  shows  the  normalized 
magnitude  of  the  kernel  of  the  integrand  (the  depth-dependent  Green's  function 
D)  of  equation  (17).  There  is  a  direct  correspondence  between  the  peaks  of 
the  kernel  in  this  figure  and  the  locations  of  the  excitation  functions  given 
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1  2  3  4  5  6  7 

123456789012345678901234567890123456789012345678901234567890123456789012345678 

&CARD1 

N0SSP=8, 

IM0DE=1, 

ITEST=0, 

IPL0T=1 , 

&END 


0.00 

1536.50 

0.00 

0.00 

0.00 

1.00 

0.00 

152.40 

1539.24 

0.00 

0.00 

0.00 

1.00 

0.00 

253.90 

1524.01 

0.00 

0.00 

0.00 

1.00 

0.00 

406.30 

1501.14 

0.00 

0.00 

0.00 

1.00 

0.00 

863.50 

1479.19 

0.00 

0.00 

0.00 

1.00 

0.00 

1015.90 

1471.88 

0.00 

0.00 

0.00 

1.00 

0.00 

5587.90 

1549.60 

0.00 

0.00 

0.00 

1.00 

0.00 

5587.90 

1555.52 

0.00 

0.00 

0.00 

1.9176 

0.00 

&CARD3 

ISRCE=3, 
IRCVR=5 , 
IPRUN=2, 
N0FFT=512, 
INEXT=0, 
ITAPR=0 , 
IDISP=0, 
IC0NT=O, 

SEND 

&CARD4 

FREQ=25. , 

SLVL=1. , 

CMIN=1471.88, 

CMAX=1555.52, 

RMIN=0.25, 

RMAX=200. , 

DELR=0.25, 

&END 

&CARD5 

REFLT0=0. , 
ALPHA0=0. , 
ATTEN0=0. , 
RH0O=O. , 

H=0. , 

SEND 


Figure  4.  Data  Input  Sequence  for  Test  Case  1 
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Sound  Speed  (m/s) 

1460  1480  1500  1  520  1  540  1  560 


Figure  5.  Sound  Speed  Profile  for  Test  Case  1  (Source  Depth  =  253.9  m, 
Receiver  Depth  =  863.5  m,  cN+1  =  1555.52  m/s,  =  1.9176  g/cm3) 
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Figure  6.  (a)  Mode  Excitation  Functions  Versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range 
Test  Case  1,  Discrete  Spectrum,  Frequency  =  25  Hz 
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FREQ  SRCE  RCVR 

17HH  MODEL  25. OO  233.39  363.50 


FREQ  S.R  C  E  RC¥R 

17HH  (FFP)  23.00  253.6*  363-30 


Figure  7.  (a)  Green's  Function  Kernel  Versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range 
Test  Case  1,  Discrete  Spectrum,  Frequency  =  25  Hz 
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in  Figure  6.  This  plot  of  the  magnitude  of  the  kernel  function,  however, 
includes  the  effect  of  the  attenuation  required  to  displace  the  poles  from  the 
real  integration  axis.  As  expected,  after  removing  the  effects  of  the  added 
attenuation,  the  propagation  loss  versus  range  curve  (figure  7(b))  for  the 
fast  field  result  agrees  almost  exactly  with  the  corresponding  curve 
determined  by  the  normal  mode  method. 

The  contribution  of  the  discrete  spectrum  gives  rise  to  those  propagating 
modes  whose  equivalent  angles  at  the  ocean  bottom  are  less  than  that  of  the 
critical  ray,  i.e.,  those  rays  near  grazing  which  reflect  from  the  ocean  floor 
without  any  loss  into  the  bottom.  This  accounts  for  the  high  losses  observed 
before  the  first  convergence  zone,  a  region  that  becomes  "filled  in"  when  the 
steeper  angle  bottom  bounce  rays  (the  continuous  spectrum)  are  included  in  the 
calculations.  This  contribution  is  accommodated  in  the  fast  field  method  by 
sampling  smaller  wavenumbers.  Figures  8  and  9  show  the  effects  of  increasing 
the  value  of  CMAX  to  include  these  wavenumbers.  In  each  case,  Ak  was  fixed  at 
the  value  used  in  figure  7  and  figure  8(a)  to  span  the  discrete  spectrum  with 
N0FFT  =  512  points.  In  figure  8(b)  the  range  of  wavenumbers  has  been  doubled 
as  N0FFT  is  increased  to  1024  points.  Figure  9  shows  the  effects  of  further 
doubling  of  the  wavenumber  range  using  values  of  N0FFT  of  2048  (figure  9(a)) 
and  4096  (figure  9(b))  points.  This  sequence  of  results  shows  the  effects  of 
increasing  steeper  angle  equivalent  rays  as  the  range  of  the  wavenumbers 
included  in  the  integration  is  increased.  Since  these  rays  correspond  to  the 
bottom  bounce  contributions,  the  region  before  the  first  convergence  zone  is 
steadily  filled  in.  While  it  is  evident  that  the  contribution  of  the  bottom 
bounce  energy  is  now  included  in  the  computation,  the  structure  of  the 
convergence  zones  are  not  affected  by  these  higher  angle  ray  arrivals.  (The 
equivalent  ray  angles  at  the  source  for  each  of  the  CMAX  values  used  in  the 
calculations  is  indicated.) 

As  a  final  computation  for  this  profile,  figure  10  shows  the  normal  mode 
result  when  the  frequency  was  increased  to  250  Hz.  In  this  case,  437 
propagating  modes  were  determined.  The  structure  of  the  excitation  functions 
versus  wavenumber  changes  significantly  from  that  shown  in  figure  5.  The 
dominant  group  of  modes  now  occurs  near  the  highest  wavenumbers.  As  a  result, 
the  propagation  loss  versus  range  curve  (figure  10(b))  shows  constructive 
interference  near  the  leading  part  of  each  zone.  NUSC17HH  can  identify  up  to 
500  propagating  modes.  This  example  was  included  to  indicate  that  apparently 
no  computational  difficulties  arise  in  determining  sound  propagation  for 
relatively  high  frequency,  deep  water  conditions. 
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FREQ  SfiCE  fiCVfi 


17HH  (FFP)  25.00  2S3..89  363.50 


FREQ 

SRCE 

1  7HH  (FFP) 

25.00 

253  av 

ES3  . 

Figure  8.  Fast  Field  Propagation  Losses  Versus  Range,  Test  Case  1 
Discrete  Plus  Continuous  Spectrum,  Frequency  =  25  Hz 

(a)  N0FFT  =  512,  CMAX  =  1555.52  m/s.  Source  Aperture  <  11.6° 

(b)  N0FFT  =  1024,  CMAX  =  1649.43  m/s.  Source  Aperture  <  22.5° 
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FREQ  SRCE  RCVR 

17HH  (FFP)  25  00  ZSJ.SS  B65  50 


FREQ  SRCE  RCVR 

17HH  (FFP)  25. QD  255.89  B5J  5D 


Figure  9.  Fast  Field  Propagation  Losses  Versus  Range,  Test  Case  1 
Discrete  Plus  Continuous  Spectrum,  Frequency  =  25  Hz 

(a)  N0FFT  =  2048,  CMAX  =  1875.95  m/s.  Source  Aperture  <  35.7° 

(b)  N0FFT  =  4096,  CMAX  =  2586.31  m/s,  Source  Aperture  <  53.9° 
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Figure  10.  (a)  Mode  Excitation  Functions  Versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range 
Test  Case  1,  Discrete  Spectrum,  Frequency  =  250  Hz 
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5.2.  TEST  CASE  2 

The  data  input  sequence  for  the  second  test  case  is  shown  in  figure  11 
and  the  sound  speed  profile  for  this  example  is  plotted  in  figure  12.  This 
profile  is  representative  of  the  Central  Arctic  Ocean  and  exhibits  a 
pronounced  surface  duct.  Since  the  entire  profile  is  upward  refracting,  the 
bottom  is  modeled  as  a  uniform  half  space  of  density  1  g/cm3  and  sound  speed 
equal  to  that  at  the  base  of  the  water  column.  For  this  example,  calculations 
were  made  at  a  frequency  of  50  Hz  for  a  source  depth  of  91.44  m  and  a  receiver 
depth  of  137.14  m. 

Using  the  fast  field  computational  mode,  the  kernel  versus  wavenumber  and 
associated  propagation  loss  versus  range  plots  are  shown  in  figure  13.  Only 
the  region  of  the  wavenumber  axis  spanned  by  the  discrete  modes  was  sampled 
with  N0FFT  =  2048  points.  The  plot  of  the  kernel  indicates  that  the 
propagation  is  dominated  by  a  single  mode  for  these  conditions.  It  is 
interesting  to  observe  that  the  wavenumber  for  this  mode  does  not  correspond 
to  an  eigenray  connecting  the  source  and  receiver,  i.e.,  the  turning  point 
depth  of  the  equivalent  ray  lies  between  the  depths  of  the  source  and  the 
receiver.  Inspection  of  the  propagation  loss  shows  that  the  character  of  the 
curve  is  well  described  by  cylindrical  spreading.  The  convergence  zone 
behaviour  of  the  deep  cycling  rays  is  apparent  at  multiples  of  about  45  km. 
Figure  14  depicts  the  kernel  and  propagation  loss  when  an  rms  surface 
roughness  of  5  m  is  included  in  the  calculations  by  setting  the  value  of 
REFLTO  =  5.  The  relative  magnitude  of  the  peaks  in  the  kernel  have  changed 
from  their  values  in  figure  13  and  the  width  of  the  dominant  peak  has 
increased.  The  effect  on  the  propagation  loss  is  to  increase  the  loss  from 
cylindrical  spreading.  This  increase  in  loss  can  be  determined  from  the  two 
graphs  to  be  approximately  0.06  dB/km. 

The  effect  of  adding  a  solid  ice  canopy  to  the  sound  speed  profile  is 
illustrated  in  figure  15.  A  5  m  thick  ice  sheet  was  inserted  at  the  top  of 
the  profile  with  a  density  of  0.9  g/cm3,  a  compressional  sound  speed  of  3500 
m/s,  and  a  shear  sound  speed  of  1800  m/s.  We  note  that  both  ISRCE  and  IRCVR 
must  be  increased  by  2  in  order  to  maintain  the  previous  source  and  receiver 
positions  within  the  water  profile.  In  addition,  all  profile  depths  were 
increased  by  5  m.  For  this  example,  since  the  wavenumbers  associated  with  the 
shear  wave  speed  were  not  expected  to  contribute  significantly  for  the  source 
and  receiver  depths  used  in  the  calculations,  only  those  wavenumbers  spanning 
the  maximum  and  minimum  values  of  compressional  sound  speed  in  the  water 
column  were  considered  in  the  integration.  As  a  result,  the  effects  of  the 
shear  waves  entered  only  through  the  boundary  conditions  satisfied  at  layer 
interfaces. 

Comparing  figures  13  and  15,  it  is  evident  that  the  presence  of  the  ice 
layer  increases  the  propagation  loss  by  only  a  few  dB.  On  the  other  hand,  the 
structure  of  the  propagation  loss  curve  in  figure  15  is  altered  considerably 
from  that  shown  in  figure  13.  Comparison  of  the  two  kernels  indicates  that 
the  ice  layer  redistributes  the  energy  between  the  lowest  three  modes.  The 
effect  on  the  propagation  loss  curve  is  to  enhance  the  interference  pattern 
between  the  second  and  third  modes  giving  rise  to  a  beating  appearance 
superimposed  on  the  cylindrical  spreading  loss  associated  with  the  dominant 
first  mode. 
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&CARD1 

N0SSP=16 , 
IM0DE=O, 
ITEST=0, 
IPL0T=1 , 

&END 


0.00 

1436.81 

0.00 

0.00 

0.00 

1.00 

0.00 

61.07 

1438.71 

0.00 

0.00 

0.00 

1.00 

0.00 

91.44 

1440.21 

0.00 

0.00 

0.00 

1.00 

0.00 

101.21 

1440.69 

0.00 

0.00 

0.00 

1.00 

0.00 

137.14 

1447.28 

0.00 

0.00 

0.00 

1.00 

0.00 

190.67 

1457.07 

0.00 

0.00 

0.00 

1.00 

0.00 

240.36 

1459.71 

0.00 

0.00 

0.00 

1.00 

0.00 

352.99 

1461.39 

0.00 

0.00 

0.00 

1.00 

0.00 

422.02 

1461.62 

0.00 

0.00 

0.00 

1.00 

0.00 

498.50 

1462.60 

0.00 

0.00 

0.00 

1.00 

0.00 

622.77 

1462.85 

0.00 

0.00 

0.00 

1.00 

0.00 

924.17 

1464.35 

0.00 

0.00 

0.00 

1.00 

0.00 

1266.52 

1467.76 

0.00 

0.00 

0.00 

1.00 

0.00 

1939.09 

1477.90 

0.00 

0.00 

0.00 

1.00 

0.00 

3657.60 

1507.14 

0.00 

0.00 

0.00 

1.00 

0.00 

3657.60 

1507.14 

0.00 

0.00 

0.00 

1.00 

0.00 

&CARD3 

ISRCE=3, 

IRCVR=5 , 

IPRUN=2, 

N0FFT=2O48, 

INEXT=0, 

ITAPR=0, 

IDISP=0, 

IC0NT=O, 

SEND 

&CARD4 

FREQ=50. , 

SLVL=1. , 

CM I N= 1436 .81 , 

CMAX=1507 . 14, 

RMIN=0.25, 

RMAX=200. , 

DELR=0.25 , 

&END 

&CARD5 

REFLT0=0. , 

ALPHA0=0. , 

ATTEN0=0. , 

RH0O=O. , 

H=0. , 

&END 

Figure  11.  Data  Input  Sequence  for  Test  Case  2. 
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Figure  13.  (a)  Green's  Function  Kernel  versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  versus  Range 
Test  Case  2,  Discrete  Spectrum,  Frequency  =  50  Hz,  a  =  0  m 
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Figure  14.  (a)  Green's  Function  Kernel  versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  versus  Range 
Test  Case  2,  Discrete  Spectrum,  Frequency  =  50  Hz,  a  =  5  m 
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Figure  15,  (a)  Greeks  Function  Kernel  versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  versus  Range 
Test  Case  2,  Discrete  Spectrum,  Frequency  =  50  Hz,  5m  Ice  Canopy 
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5.3.  TEST  CASE  3 

The  data  input  sequence  for  test  case  3  is  given  in  figure  16.  For  this 
shallow  water  profile,  shown  in  figure  17,  the  bottom  is  expected  to  play  an 
important  role  in  the  propagation  of  underwater  sound.  Moreover,  the  effects 
of  the  continuous  spectrum  will  be  seen  to  be  important.  Calculations  were 
carried  out  for  a  frequency  of  100  Hz,  a  source  depth  of  30  m,  and  a  receiver 
depth  of  90  m. 

Figure  18  shows  the  discrete  mode  results  for  this  test  case.  Only  three 
proper  modes  were  able  to  propagate  under  these  conditions.  The  coherent  sum 
of  the  three  modes  (figure  18(a))  determines  the  propagation  loss  versus  range 
curve  shown  in  the  bottom  part  of  the  figure.  With  only  three  modes 
propagating,  the  character  of  the  propagation  loss  versus  range  curve  (figure 
18(b))  is  seen  to  be  smooth.  Figure  19  shows  the  results  when  the  fast  field 
computational  option  is  used  and  part  of  the  continuous  spectrum  is  included 
in  the  integration  (CMAX  =  1600  m/s).  It  is  apparent  from  the  plot  of  the 
kernel  (figure  19(a))  that  in  addition  to  the  three  peaks  at  the  locations  of 
the  wavenumbers  of  the  discrete  modes,  a  considerable  contribution  arises  in 
the  adjacent  span  of  the  continuous  spectrum.  The  effect  on  the  propagation 
loss  curve  (figure  19(b))  is  to  modulate  the  discrete  mode  result  of  the 
previous  figure. 

According  to  the  flowchart  of  figure  3,  an  alternate  option  of  obtaining 
the  fast  field  result  is  to  evaluate  numerically  the  branch  line  integral 
associated  with  the  modal  representation.  Figure  20  presents  the  result  of 
this  computation  when  the  branch  line  integral  is  evaluated  via  the  fast  field 
method  by  setting  the  flag  IC0NT  <  0.  N0FFT  was  set  to  4096  points  to  ensure 
that  the  span  of  the  unaliased  range  values  reached  20  km.  For  comparison, 
another  fast  field  computation  is  shown  in  figure  20(a)  where  CMAX  was 
increased  to  1800  m/s  and  N0FFT  was  increased  to  2048  points.  While  the 
propagation  loss  versus  range  curves  are  similar  for  the  two  methods  of 
computation,  the  differences  can  be  attributed  to  the  points  raised  in  the 
discussion  in  section  3.5.  In  particular,  the  branch  line  integral  only 
included  the  portion  around  the  real  axis  and  the  sound  speed  profile  of  the 
water  column  was  taken  to  consist  of  isovelocity  layers  for  the  branch  line 
evaluation.  It  is  apparent  that  the  fast  field  method  of  calculation  gives 
the  correct  resuU  with  much  less  effort  in  this  case.  Efforts  to  compute  the 
branch  line  contribution  via  numerical  integration  using  Gaussian  quadrature 
techniques  did  not  prove  successful. 
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&CARD1 

N0SSP=6 , 

IM0DE=1 , 

ITEST=0, 

I PL0T= 1 , 

&END 


0.00 

1500.00 

0.00 

0.00 

0.00 

1.00 

0.00 

30.00 

1499.50 

0.00 

0.00 

0.00 

1.00 

0.00 

90.00 

1498.50 

0.00 

0.00 

0.00 

1.00 

0.00 

120.00 

1498.00 

0.00 

0.00 

0.00 

1.00 

0.00 

240.00 

1500.00 

0.00 

0.00 

0.00 

1.00 

0.00 

240.00 

1505.00 

0.00 

0.00 

0.00 

2.10 

0.00 

&CARD3 

ISRCE=2, 
IRCVR=3, 
IPRUN=2, 
N0FFT=1O24, 
INEXT=0, 
ITAPR=0, 
IDISP=0 , 
IC0NT=O , 

&END 

&CARD4 

FREQ=100. , 
SLVl=l. , 
CMIN=1498. , 
CMAX=1505. , 
RMIN=0. 1 , 
RMAX=20. , 
DELR=0. 1 , 

&END 

&CARD5 

REFLT0=0. , 
ALPHA0=0. , 
ATTEN0=0. , 
RH0O=O. , 
H=0. , 

&END 


Figure  16.  Data  Input  Sequence  for  Test  Case  3. 
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FREQ  SRCE  RCVR 


1  7HH  MODEL  1  00  00  30.00  90  00 


‘390  395  400  405  41  0  41  5  420  425 


Wavenumber  (/km) 

FREQ  SRCE  RCVR 

17HH  (MODE)  100.00  30  00  90.00 


Range  (km) 


Figure  18.  (a)  Mode  Excitation  Functions  versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range 
Test  Case  3,  Discrete  Spectrum,  Frequency  =  100  Hz 
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FREQ  SRCE  RCVR 


17HH  MODEL  100.00  30.00  90  00 


FREQ  SR Cl  RCVR 


17HH  (FFP)  100.00  30.00  90 . 00 


Range  (km) 


Figure  19.  (a)  Green's  Function  Kernel  Versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range 
Test  Case  3,  N0FFT  =  1024,  CMAX  =  1650  m/s.  Frequency  =  100  Hz 


37 


Ltiss  (  d  B 


TR  7221 


ffteo 

I  00 .  DO 


SRCE 
30  .00 


RCVR 
90 .00 


17HH  (FFP) 


FREO  SRCE  RCVR 


1  7  HH  (MODE)  100  00  30  00  90.00 


Figure  20.  Propagation  Losses  versus  Range,  Test  Case  3,  Frequency  =  100  Hz 
(a)  Fast  Field  Method,  N0FFT  =  2048,  CMAX  =  1800  m/s 
(b)  Normal  Mode  Method,  N0FFT  =  4096,  IC0NT  <  0 
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5.4.  TEST  CASE  4 

This  final  example  was  taken  from  a  recent  workshop17  and  was  included 
here  to  represent  propagation  in  shallow  water  over  a  lossy  bottom.  An 
excellent  account  of  all  the  test  cases  in  the  workshop  has  been  carried  out 
using  another  version  of  the  fast  field  model.18  The  data  input  sequence  for 
this  example  is  given  in  figure  21.  This  example  is  essentially  an  extension 
of  the  so-called  Pekeris  model  in  which  the  absorption  of  the  bottom  is  not 
zero.  The  sound  speed  profile  in  figure  22  shows  an  isovelocity  water  layer 
of  sound  speed  1500  m/s  and  density  1  g/cm3  overlaying  a  lossy  fluid  halfspace 
of  sound  speed  1590  m/s,  density  2.1  g/cm3  and  absorption  76.83  dB/km  (0.5 
dB/wavel ength) .  The  calculations  were  carried  out  for  a  frequency  of  250  Hz. 

Since  there  is  attenuation  in  the  bottom  half-space  for  this  example,  the 
normal  mode  method  of  computation  cannot  be  used.  Figure  23  shows  the  results 
of  a  fast  field  evaluation  with  the  source  and  receiver  located  in  the  middle 
of  the  water  layer  at  a  depth  of  50  m.  Figure  23(a)  presents  the  magnitude  of 
the  kernel  as  a  function  of  wavenumber  in  which  11  propagating  modes  can  be 
identified.  For  this  source  and  receiver  depth,  it  is  evident  that  the  lower 
order  modes  (slowest  phase  velocities)  dominate  the  propagation.  The  effect 
of  attenuation  is  not  readily  apparent  on  the  propagation  loss  versus  range 
curve  (figure  23(b))  since  the  lower  order  modes  do  not  penetrate 
significantly  into  the  bottom  half-space.  Only  the  range  interval  between  5 
km  and  10  km  has  been  plotted.  Within  this  range  interval  ,  only  the 
propagating  modes  contribute  to  the  acoustic  field,  thus  only  that  region  of 
the  wavenumber  axis  was  sampled  in  the  integration. 

The  effect  of  locating  the  source  and  receiver  near  the  bottom  of  the 
water  layer  at  a  depth  of  99.5  m  is  shown  in  figure  24.  The  distribution  of 
dominant  modes  has  shifted  to  the  higher  mode  numbers  (figure  24(a)).  As  a 
result,  the  acoustic  field  within  the  water  is  influenced  more  by  the 
properties  of  the  lower  half-space  since  the  higher  order  modes  penetrate 
further  into  the  bottom.  The  large  increase  in  propagation  loss  (figure 
24(b))  observed  near  a  range  of  7  km  is  due  to  the  presence  of  absorption  in 
the  bottom  and  its  considerable  effect  on  these  more  deeply  penetrating  higher 
order  modes. 
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1  2  3  4  5  6  7 

123456789012345678901234567890123456789012345678901234567890123456789012345678 

&CARD1 

N0SSP=4, 

IM0DE=O, 

ITEST=0, 

IPL0T=1, 

&END 


0.00 

1500.00 

0.00 

0.00 

0.00 

1.00 

0.00 

50.00 

1500.00 

0.00 

0.00 

0.00 

1.00 

0.00 

100.00 

1500.00 

0.00 

0.00 

0.00 

1.00 

0.00 

100.00 

1590.00 

0.00 

76.83 

0.00 

1.20 

0.00 

&CARD3 

ISRCE=2, 

IRCVR=2, 

IPRUN=2, 

N0FFT=512 , 

INEXT=0, 

ITAPR=0, 

IDISP=0, 

IC0NT=O, 

SEND 

&CARD4 

FREQ=250. , 
SLVL=1. , 
CMIN=1490. , 
CMAX=1600. , 
RMIN=0.1, 
RMAX=10. , 
DELR=0.1, 

&END 

&CARD5 

REFLT0=0. , 
ALPHA0=0. , 
ATTEN0=0. , 
RH0O=O. , 
H=0. , 

SEND 


Figure  21.  Data  Input  Sequence  for  Test  Case  4. 
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Sound  Speed  (m/s) 


o_ 

0 

o 


1  450 

o 


20  - 


40 


60 


80 


1  00 


1  20 


1  500 

—Q — 


1  550 


1  600 


T 


Figure  22.  Sound  Speed  Profile  for  Test  Case  4 
(c^+2  =  1590  m/s,  =  1.2  g/cin^,  aN+l.  =  dB/km) 
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FREQ  SRCE  RCVR 

17HH  MODEL  250.00  50.00  50.00 


FREQ  SRCE  RCVR 

1 7HH  (FFP)  250.00  50.00  50.00 
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Figure  23.  (a)  Green's  Function  Kernel  Versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range,  Test  Case  4 
Source  Depth  =  Receiver  Depth  =  50.0  m,  Frequency  =  250  Hz 
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FREQ  SRCE  RCVR 


1  7  H  H  MODEL  250,00  99  50  99.50 


FREQ  SRCE  RCVR 

17HH  (FFP)  250  00  99  50  99  50 


Figure  24.  (a)  Green's  Function  Kernel  Versus  Wavenumber,  and 

(b)  Associated  Propagation  Losses  Versus  Range,  Test  Case  4 
Source  Depth  =  Receiver  Depth  =  99.5  m.  Frequency  =  250  Hz 
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6.  SUMMARY 


This  report  describes  the  implementation  of  the  17HH  version  of  the  LDGO 
normal  mode/fast  field  model  on  the  NUSC  VAX  780/11  computer.  This  model 
differs  from  most  propagation  models  used  in  underwater  acoustics  since  it  can 
treat  both  compressional  and  shear  wave  propagation  within  a  layered 
fluid/solid  medium.  Sufficient  documentation  is  provided  to  enable 
prospective  NUSC  users  to  obtain  predictions  of  sound  propagation.  A  plotting 
program  is  also  described  for  obtaining  quality  graphs  generated  by  the  normal 
mode  or  fast  field  methods  of  computation.  Several  test  cases  have  been 
provided  (complete  with  data  input  sequences)  to  illustrate  the  use  of  the 
model  under  representative  environmental  conditions.  Further  examples  of  fast 
field  method  calculations  may  be  found  in  a  report  by  Jensen  and  Kuperman18  in 
which  several  propagation  loss  models  are  compared. 

The  execution  times  for  each  of  the  test  cases  vary  from  a  few  seconds  to 
more  than  an  hour.  Table  1  summarizes  the  actual  CPU  times  required  to  obtain 
the  information  displayed  in  each  Figure. 


Table  1.  CPU  Execution  Times  for  the  Test  Cases. 


Test  Case 

Fi gure 

CPU  (minutes) 

1 

6 

8.7 

7 

3.5 

8(a) 

3.5 

8(b) 

8.1 

9(a) 

16.5 

9(b) 

31.5 

10 

63.9 

2 

13 

17.9 

14 

17.1 

15 

18.9 

3 

18 

3.1 

19 

5.5 

20(a) 

10.1 

20(b) 

12.5 

4 


23 

24 


0.5 

0.5 
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APPENDIX 


EQUIVALENCE  LIST  BETWEEN  LDGO  AND  NUSC  DATA  INPUT  SEQUENCES 


In  the  NUSC  version  of  17HH,  the  data  input  sequence  has  been  modified 
from  that  supplied  by  LDGO.  In  addition  to  making  use  of  the  FORTRAN  NAMELIST 
feature,  most  of  the  input  variables  have  been  renamed  and  their  order  of 
appearance  within  the  data  file  has  been  changed.  A  comparison  between  the 
two  input  lists  is  given  below.  Upper  case  letters  denote  variable  names 
which  appear  as  indicated  on  each  card  image.  Lower  case  letters  identify 
equivalent  variable  names  which  appear  on  different  card  images.  For  example, 
KSF  in  the  LDGO  version  has  been  renamed  ISRCE  in  the  NUSC  version  and  both 
appear  on  the  third  card  image.  On  the  other  hand,  PO  in  the  LDGO  version  has 
been  renamed  SLVL  but  has  been  relocated  from  the  third  card  image  to  the 
fourth  card  image. 

LDGO  VERSION  NUSC  VERSION 


CARD  1  F0RMAT(2I3)  NAMELIST  &CARD1 


MMAXF 

NCARD 

Not  available 
Not  available 


N0SSP 

IPL0T 

IM0DE 

ITEST 


CARD  2  ((V(I,J),J=1,7), 1=1 ,MMAXF) 


CARD2  ((V(I,J),J=1,7) ,  1  =  1  ,N0SSP) 
F0RMAT(7D1O.O) 

ct  &  3  in  dB/km 


F0RMAT(7D1O.O) 
a  &  3  in  Np/m 


CARD  3  F0RMAT(7I5 ,4D10.0) 


NAMELIST  &CARD3 


KSF 

LRF 

MM 

LS 

KKK 


ISRCE 

IRCVR 

IPRUN 

N0FFT 

INEXT 

ITAPR 


IBEAM 


IMAX 

FREQ 

PO 

C 

CMAX 


(rmax-rmin)/del r  see  &CARD4 


freq 

slvl 

cmin 

cmax 


see  &CARD4 
see  &CARD4 
see  &CARD4 
see  &CARD4 


khpv  see  CARD4 

ibr  see  CARD4 


IDISP 

IC0NT 
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LDGO  version  NUSC  version 


F0RMAT(5D10. 0,515) 

NAMELIST  &CARD4 

REFLTO 

ref I tO 

see 

&CARD5 

ALPHAO 

al phaO 

see 

&CARD5 

ATTENO 

attenO 

see 

&CARD5 

RH0O 

rhoO 

see 

&CARD5 

H 

h 

see 

&CARD5 

INAUT 

Set  to 

0 

internal ly 

L0NG 

Set  to 

1 

internal ly 

KHPV 

idi  sp 

see 

&CARD3 

IBR 

icont 

see 

&CARD3 

IDPL 

Set  to 

1 

internal ly 

freq 

see  CARD3 

FREQ 

pO 

see  CARD3 

SLVL 

c 

see  CARD3 

CMIN 

cmax 

see  CARD3 

CMAX 

ra 

see  CARD5 

RMIN 

ra+imax*di ra 

RMAX 

dl  ra 

see  CARD5 

DELR 

CARD  5  F0RMAT(2D1O.O)  NAMELIST  &CARD5 


RA 

rmin 

see  &CARD4 

DLRA 

del  r 

see  &CARD4 

refltO 

see  CARD4 

REFLTO 

al phaO 

see  CARD4 

ALPHAO 

attenO 

see  CARD4 

ATTENO 

rhoO 

see  CARD4 

RH0O 

h 

see  CARD4 

H 
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